A marine heatwave in the Santa Barbara Channel: impacts on kelp forest canopies and ecosystem state
Sofia Urgoiti Crespo
2023-10-05
Background
Marine heatwaves (MHWs) are usually defined as a coherent area of extreme warm sea surface temperature (SST) persisting for days to months. MHWs are becoming more frequent and intense, a trend which is likely to accelerate with global warming (Frolicher et al., 2018). The scientific community is still seeking to understand the processes leading to their build-up, persistence and decay. While some ecological impacts of marine heatwaves are perceived more easily, other indirect effects might affect ecosystems in the longer term. Warmer water is poorer nutrients, holds less oxygen, and has been associated with increased toxic algal blooms. Another direct impact of marine heatwaves can be massive bird die-offs due to a lack of forage fish.
Perhaps a major cause of rippling effects due to increased temperatures is a reduced kelp forest canopy. Kelp can only thrive in temperate climates where temperature do not exceed 22-24°C (Becheler et al., 2022), and exposure to heat stress has been thought to erode its capacity for canopy regeneration via reduced recruitment (Wernberg et al., 2010). Kelp forests are keystone species: they create habitat for thousands of species, and provide a wide array of ecosystem services to society. Loss of kelp could lead to reduced habitat, food availability, changes in light exposure and growth of other algae. In some cases, if purple urchin populations are widespread, they can overgraze and suppress kelp recovery, sometimes forming areas completely deprived of kelp, also known as urchin barrens (Bennett & Catton, 2019).
In late October 2013, a sea surface temperature anomaly was detected in the northeast Pacific Ocean, soon becoming the largest MHW (often called “The Blob”) ever recorded. During “The Blob”, lasting between 2013 and 2015 maximum SST anomalies of >6 °C were recorded off the coast of southern California. The impacts of “the Blob” soon became a major focus of interest among scientists, yet up until today many questions remain unanswered. A recent paper by Daniel Reed and other scientists, confirmed the resilience of certain kelp forests in mainland Santa Barbara (Reed et al., 2016), while assessments of the impacts on the Channel Islands remain scarce.
Impact does not only mean canopy loss: it can also translate into altered functional relationships between organisms or even changes in ecosystem state. This is an aspect that I aim to explore in my current scientific research as a Master’s student at UCSB. I want to investigate the following questions: (1) how resilient have Santa Barbara’s kelp forests been agains disturbance, and what is their population recovery pathway? (2) Are there any biogeographic differences in population trajectory? (3) Most importantly, I want to record any changes in ecosystem state (forested, mixed, urchin barrens) and develop a model that is able to predict
This coding exercise is part of a broader research project, which investigates the resilience of Santa Barbara’s (CA) kelp forests to temperature disturbances, as well as the influence of different drivers on kelp population recovery. Here, I will produce some maps to display the region of interest in my thesis, as well as visually explore some of the variables that I am interested in (SST and bathymetry depth) as physical drivers of the system.
Libraries
library(here)## Warning: package 'here' was built under R version 4.3.1## here() starts at C:/Users/sofiaurgoiti/Documents/R-PROJECTS-UCSB/EDS223-GEOSPATIAL-ANALYSIS/Marine-Heatwave-Maps-Final-Portfoliolibrary(devtools)## Warning: package 'devtools' was built under R version 4.3.2## Loading required package: usethis## Warning: package 'usethis' was built under R version 4.3.2library(sf)## Warning: package 'sf' was built under R version 4.3.1## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUElibrary(raster)## Warning: package 'raster' was built under R version 4.3.1## Loading required package: sp## Warning: package 'sp' was built under R version 4.3.1## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.library(stars)## Warning: package 'stars' was built under R version 4.3.1## Loading required package: abindlibrary(terra)## Warning: package 'terra' was built under R version 4.3.1## terra 1.7.55library(tmap)## Warning: package 'tmap' was built under R version 4.3.2## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.1##
## Attaching package: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(stringr)
library(geomtextpath)## Warning: package 'geomtextpath' was built under R version 4.3.2## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.3.1library(ggplot2)
library(ggspatial)## Warning: package 'ggspatial' was built under R version 4.3.1library(gt)## Warning: package 'gt' was built under R version 4.3.1library(RColorBrewer)
library(viridis)## Warning: package 'viridis' was built under R version 4.3.1## Loading required package: viridisLiteLoading the data:
The title of the dataset available from the NOAA-ERDDAP data server is “SST, Aqua MODIS, NPP, 0.0125°, West US, Day time (11 microns), 2002-present (Monthly Composite)”. The data was accessed via this form (https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMWsstdmday.html), and the SST variable used is a monthly composite for the period of 2010-2022. Data was extracted at approximately 1.47km resolution, and SST measurements accurate to within ±1 degree Celsius.
here()## [1] "C:/Users/sofiaurgoiti/Documents/R-PROJECTS-UCSB/EDS223-GEOSPATIAL-ANALYSIS/Marine-Heatwave-Maps-Final-Portfolio"# Site Location
sites_data <- read.csv(here("Data", "upc_swath_ucsb.csv"))
#SST
#pre heatwave: 2010,2011,2012, 2013
sst_2010_06 <- rast(here("Data", "sst_2010_06.tif"))
sst_2011_06 <- rast(here("Data", "sst_2011_06.tif"))
sst_2012_06 <- rast(here("Data", "sst_2012_06.tif"))
sst_2013_06 <- rast(here("Data", "sst_2013_06.tif"))
#during heatwave: 2014, 2015, 2016
sst_2014_06 <- rast(here("Data", "sst_2014_06.tif"))
sst_2015_06 <- rast(here("Data", "sst_2015_06.tif"))
sst_2016_06 <- rast(here("Data", "sst_2016_06.tif"))
#post heatwave: 2017, 2018, 2019, 2020, 2021
sst_2017_06 <- rast(here("Data", "sst_2017_06.tif"))
sst_2018_06 <- rast(here("Data", "sst_2018_06.tif"))
sst_2019_06 <- rast(here("Data", "sst_2019_06.tif"))
sst_2020_06 <- rast(here("Data", "sst_2020_06.tif"))
sst_2021_06 <- rast(here("Data", "sst_2021_06.tif"))
sst_2022_06 <- rast(here("Data", "sst_2022_06.tif"))
# Why are the SST values unusually high? They do not make sense ...
summary(sst_2010_06)## Warning: [summary] used a sample## sst_2010_06
## Min. : 0.00
## 1st Qu.: 0.00
## Median : 86.00
## Mean : 80.33
## 3rd Qu.:138.00
## Max. :241.00summary(sst_2010_06, maxpixels = ncell(sst_2010_06))## Warning: [summary] used a sample## sst_2010_06
## Min. : 0.00
## 1st Qu.: 0.00
## Median : 86.00
## Mean : 80.33
## 3rd Qu.:138.00
## Max. :241.00min(values(sst_2010_06))## [1] 0max(values(sst_2010_06))## [1] 255mean(values(sst_2010_06))## [1] 80.27975# Bathymetry
bathymetry_WC <- rast(here("Data", "depth.tif")) #data from GEBCO Gridded Bathymetry data
bathymetry_SB <- rast(here("Data", "gebco_2022_n34.7291_s33.5033_w-120.9447_e-119.048.tif")) #GEBCO gridded bathymetry data, cropped at the SBC scaleWrangling data:
# Subsetting site coordinates in dataframe
sites_coords <- sites_data %>%
select(campus, survey_year, site, site_status, longitude, latitude) %>%
relocate(site, .before = longitude)
# Transform site coordinates to an sf_object with a coordinate system so we can plot them
sites_sf <- st_as_sf(sites_coords, coords = c("longitude", "latitude"), crs = 4326) %>%
relocate(site, .before = "campus")
# Stacking SST data
sst_stack <- c(sst_2010_06, sst_2011_06, sst_2012_06, sst_2013_06, sst_2014_06, sst_2015_06, sst_2016_06, sst_2017_06, sst_2018_06, sst_2019_06, sst_2020_06, sst_2021_06, sst_2022_06) #using terra package format to stack rasters: the stacking worked, which indicates that rasters are of the same extent, resolution, and CRS
#Checking CRS coordinates of datasets & check they loaded correctly
crs(bathymetry_SB)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"crs(sst_stack)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"crs(sites_sf)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"# Visually check that they loaded correctly
plot(bathymetry_WC)
plot(bathymetry_SB)
plot(sst_2010_06) #here we only look at one layer
plot(sst_stack)
1. Defining the spatial context: Santa Barbara Channel and PISCO sites
First I want to make some simple maps of both regions of interest: the mainland Santa Barbara coastline and the Channel Islands.
# Map of sites only: ---------------------------------------------------------
# Define the bounding box (xmin, xmax, ymin, ymax)
bbox_isl <- c(xmin = -120.5, xmax = -119.25, ymin = 33.8, ymax = 34.10) #bbox for the islands
bbox_main <- c(xmin = -120.5, xmax = -119.5, ymin = 34.3, ymax = 34.5)
# Create an EXTENT object for each
extent_isl <- ext(bbox_isl[1], bbox_isl[2], bbox_isl[3], bbox_isl[4])
extent_main <- ext(bbox_main[1], bbox_main[2], bbox_main[3], bbox_main[4])
# Crop the raster using the bounding box
zoomed_island <- crop(bathymetry_SB, extent_isl)
zoomed_mainland <- crop(bathymetry_SB, extent_main)
# Plot the original and cropped raster
par(mfrow = c(1, 2))
plot(zoomed_island, main = "Cropped Islands Raster")
plot(zoomed_mainland, main = "Cropped Mainland Raster")
#Plotting both:
tmap_mode("plot")## tmap mode set to plottingtm_shape(zoomed_island) +
tm_raster() +
tm_shape(sites_sf) +
tm_squares(col = "black", size = 0.05) ## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
2. Bathymetry Data
#1. Basic GGPLOT ----------------------------------------
#convert the raster to a dataframe. This means we extract the x and y coordinates from the raster data.
bathymetry_df <- data.frame(values = values(bathymetry_SB), xyFromCell(bathymetry_SB, 1:ncell(bathymetry_SB)))
#rename column names
colnames(bathymetry_df) <- c("layer", "lon", "lat")
# Create a color palette for depth
blue_palette <- colorRampPalette(c("lightcyan", "royalblue4"))
shades_blue <- blue_palette(5)# Generate 8 shades of blue# Bathymetry Map 1 ---------------------------------------
ggplot(bathymetry_df, aes(x = lon, y = lat, fill = layer)) +
geom_raster() +
scale_fill_gradientn(colors = rev(shades_blue), limits = c(-5000, 0)) +
theme_minimal() +
theme(
panel.grid = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
) +
labs(title = "Santa Barbara Channel Bathymetry Map", x = "", y = "", fill = "Depth (m)") +
annotation_scale(location = "br", scale = 100000, pad_x = unit(0.4, "in"), pad_y = unit(0.3, "in")) + #scale parameter is set to 100,000, since our data is in m and we want the scale bar to display km
annotation_north_arrow(location = "br", pad_x = unit(0.3, "in"), pad_y = unit(0.5, "in"), style = ggspatial::north_arrow_nautical(fill = c("black", "white")))## Warning in annotation_scale(location = "br", scale = 1e+05, pad_x = unit(0.4, :
## Ignoring unknown parameters: `scale`## Using plotunit = 'm'
#2. Bathymetry Map 2: Adding depth contour lines ------------------------------------------
bathymetry_SB[bathymetry_SB > 0] <- NA #set depths above 0 to NA so they don't appear on the map
contour_lines <- contour(bathymetry_SB) #generate lines
contour_df <- data.frame(xyFromCell(bathymetry_SB, 1:ncell(bathymetry_SB)), depth = values(bathymetry_SB)) #extract values and coordinates
colnames(contour_df) <- c("lon", "lat", "depth") #renaming columns
#Plot bathymetry (in blue) & contours
ggplot() +
geom_raster(data = bathymetry_df, aes(x = lon, y = lat, fill = layer)) +
scale_fill_gradientn(colors = rev(shades_blue), limits = c(-5000, 0)) +
geom_contour(data = contour_df, aes(x = lon, y = lat, z = depth), color = "black", n = 1, size = 0.2, alpha = 0.5) +
theme_minimal() +
theme(plot.title = element_text(size = 8),
panel.grid = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
) +
labs(title = "Santa Barbara Channel Bathymetry Map", x = "", y = "", fill = "Depth (m)") +
annotation_scale(location = "tr", line_width = 0.5, height = unit(0.1, "cm"), text_cex = 0.5, scale = 100000, pad_x = unit(0.35, "in"), pad_y = unit(0.2, "in")) + #scale parameter is set to 100,000, since our data is in m and we want the scale bar to display km
annotation_north_arrow(location = "tr", height = unit(1, "cm"), width = unit(1, "cm"), pad_x = unit(0.2, "in"), pad_y = unit(0.3, "in"), style = ggspatial::north_arrow_nautical(fill = c("black", "white")))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning in geom_contour(data = contour_df, aes(x = lon, y = lat, z = depth), :
## Ignoring unknown parameters: `n`## Warning in annotation_scale(location = "tr", line_width = 0.5, height =
## unit(0.1, : Ignoring unknown parameters: `scale`## Warning: Removed 34357 rows containing non-finite values (`stat_contour()`).## Using plotunit = 'm'
# trying to add text, but I cannot find a way to make it less crowded
ggplot() +
geom_raster(data = bathymetry_df, aes(x = lon, y = lat, fill = layer)) +
scale_fill_gradientn(colors = rev(shades_blue), limits = c(-5000, 0)) +
geom_contour(data = contour_df, aes(x = lon, y = lat, z = depth), color = "black", n = 1, size = 0.2, alpha = 0.5) +
geom_text(data = subset(contour_df, depth %% 500 == 0), #trying to add labels next to contour lines
aes(x = lon, y = lat,
label = round(depth, digits = 1)),
nudge_x = 0.03, nudge_y = 0.01, #controls location of labels
size = 2, color = "black") +
theme_minimal() +
theme(plot.title = element_text(size = 8),
panel.grid = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
) +
labs(title = "Santa Barbara Channel Bathymetry Map", x = "", y = "", fill = "Depth (m)") +
annotation_scale(location = "tr", line_width = 0.5, height = unit(0.1, "cm"), text_cex = 0.5, scale = 100000, pad_x = unit(0.35, "in"), pad_y = unit(0.2, "in")) + #scale parameter is set to 100,000, since our data is in m and we want the scale bar to display km
annotation_north_arrow(location = "tr", height = unit(1, "cm"), width = unit(1, "cm"), pad_x = unit(0.2, "in"), pad_y = unit(0.3, "in"), style = ggspatial::north_arrow_nautical(fill = c("black", "white")))## Warning in geom_contour(data = contour_df, aes(x = lon, y = lat, z = depth), :
## Ignoring unknown parameters: `n`## Warning in annotation_scale(location = "tr", line_width = 0.5, height =
## unit(0.1, : Ignoring unknown parameters: `scale`## Warning: Removed 34357 rows containing non-finite values (`stat_contour()`).## Using plotunit = 'm'
#3. Interactive map version -------------------------------------
tmap_mode("view")## tmap mode set to interactive viewingbathy_interactive <- tm_shape(bathymetry_SB) +
tm_raster(palette = shades_blue) +
tm_shape(sites_sf) +
tm_bubbles(col = "black", size = 0.025) +
tm_layout(legend.outside = TRUE,
main.title.size = 1,
main.title = "Santa Barbara Channel Bathymetry Map",
frame = T) +
tm_compass(type = "arrow",
position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "top"))
bathy_interactive## Compass not supported in view mode.# Now make the same map but even more zoomed in the islands/mainland
tmap_mode("plot")## tmap mode set to plottingtm_shape(zoomed_mainland) +
tm_raster(palette = shades_blue) +
tm_shape(sites_sf) +
tm_bubbles(col = "black", size = 0.05) +
tm_layout(legend.outside = TRUE,
main.title.size = 1,
main.title = "Santa Barbara Channel Bathymetry Map",
frame = T) +
tm_compass(type = "arrow",
position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "bottom"))
#
tm_shape(zoomed_island) +
tm_raster(palette = shades_blue) +
tm_shape(sites_sf) +
tm_bubbles(col = "black", size = 0.05) +
tm_layout(legend.outside = TRUE,
main.title.size = 1,
main.title = "Santa Barbara Channel Bathymetry Map",
frame = T) +
tm_compass(type = "arrow",
position = c("left", "bottom")) +
tm_scale_bar(position = c("right", "bottom"))
3. Surface Temperature Data
#1.
# We now want to crop the SST raster stack to our region of interest (from the whole of the WC to the SBC) - to do so, we crop the sst raster with the bathymetry_SB
sst_stack <- project(sst_stack, bathymetry_SB) #set the CRS and extent of sst to match bathymetry raster
# Crop SST raster to match the extent of bathymetry raster
sst_stack_crop <- crop(sst_stack, bathymetry_SB)
#Values Check
summary(sst_stack_crop) #checking range of temperature values in all layers## Warning: [summary] used a sample## sst_2010_06 sst_2011_06 sst_2012_06 sst_2013_06
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 29.26 1st Qu.: 13.17 1st Qu.: 14.17 1st Qu.: 15.56
## Median : 86.08 Median : 92.33 Median : 95.00 Median : 93.11
## Mean : 71.56 Mean : 73.36 Mean : 78.45 Mean : 78.39
## 3rd Qu.:104.00 3rd Qu.:104.75 3rd Qu.:114.25 3rd Qu.:114.00
## Max. :125.67 Max. :125.56 Max. :137.67 Max. :141.03
## sst_2014_06 sst_2015_06 sst_2016_06 sst_2017_06
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 19.72 1st Qu.: 12.99 1st Qu.: 15.47 1st Qu.: 0.00
## Median :124.00 Median :106.97 Median : 81.92 Median : 88.33
## Mean : 98.35 Mean : 84.44 Mean : 66.11 Mean : 73.79
## 3rd Qu.:139.25 3rd Qu.:118.72 3rd Qu.: 95.03 3rd Qu.:107.75
## Max. :161.53 Max. :145.22 Max. :121.42 Max. :149.81
## sst_2018_06 sst_2019_06 sst_2020_06 sst_2021_06
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 12.65 1st Qu.: 15.59
## Median : 79.86 Median : 0.00 Median : 72.67 Median : 80.06
## Mean : 68.19 Mean : 45.05 Mean : 60.94 Mean : 67.06
## 3rd Qu.:100.78 3rd Qu.: 93.58 3rd Qu.: 88.92 3rd Qu.: 97.92
## Max. :133.17 Max. :146.11 Max. :111.69 Max. :117.67
## sst_2022_06
## Min. : 0.00
## 1st Qu.: 16.11
## Median : 75.31
## Mean : 63.80
## 3rd Qu.: 92.50
## Max. :113.00# Visual Check
plot(sst_stack_crop) #now we can see all the layers
plot(sst_stack_crop, col = rev(terrain.colors(10)))
plot(sst_stack_crop[[2]], col = rev(terrain.colors(10))) #select an individual layer, in this case the year 2011
Below, I develop a function that creates a map for each year from 2010-2022
#Create raster list:
raster_names <- c("sst_2010_06", "sst_2011_06", "sst_2012_06", "sst_2013_06", "sst_2014_06", "sst_2015_06", "sst_2016_06", "sst_2017_06", "sst_2018_06", "sst_2019_06", "sst_2020_06", "sst_2021_06", "sst_2022_06")
# Create a list with raster objects and their names
sst_raster_list <- mget(raster_names)
# Create the function to create a map for each year of the raster
create_sst_maps <- function(raster_list, crop_raster)
{
for (i in seq_along(raster_list))
{current_raster <- raster_list[[i]]
#we crop it to the SBC region
current_raster <- crop(current_raster, crop_raster)
#now we turn this raster into a dataframe
sst_df <- data.frame(values = values(current_raster), xyFromCell(current_raster, 1:ncell(current_raster)))
colnames(sst_df) <- c("sst", "lon", "lat") }
#now we make plot for the one layer (one year SST data)
plot <- ggplot(sst_df, aes(x = lon, y = lat, fill = sst)) +
geom_raster() +
scale_fill_gradientn(colors = terrain.colors(10)) +
theme_minimal() +
theme(
panel.grid = element_blank(), # Remove major grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
) +
labs(title = paste("Santa Barbara SST Map", i), x = "", y = "", fill = "SST (Degrees Celsius)") +
annotation_scale(location = "br", scale = 100000, pad_x = unit(0.4, "in"), pad_y = unit(0.3, "in")) + #scale parameter is set to 100,000, since our data is in m and we want the scale bar to display km
annotation_north_arrow(location = "br", pad_x = unit(0.3, "in"), pad_y = unit(0.5, "in"), style = ggspatial::north_arrow_nautical(fill = c("black", "white")))
#now we print each of the plots
print(plot)
}
create_sst_maps(sst_raster_list, bathymetry_SB)## Warning in annotation_scale(location = "br", scale = 1e+05, pad_x = unit(0.4, :
## Ignoring unknown parameters: `scale`## Using plotunit = 'm'
References
“The Blob” National Park Service https://www.nps.gov/articles/theblob.htm#:~:text=The%20Blob%20Appears&text=This%20mass%20of%20water%20was,masses%20between%202013%20and%202018
Jacox, M. G. et al. Impacts of the 2015–2016 El Nino on the California Current System: early assessment and comparison to past events. Geophys. Res. Lett. 43, 7072–7080 (2016).
Frölicher, T.L., Fischer, E.M. and Gruber, N., 2018. Marine heatwaves under global warming. Nature, 560(7718), pp.360-364.
Wernberg, T., Thomsen, M.S., Tuya, F., Kendrick, G.A., Staehr, P.A. and Toohey, B.D., 2010. Decreasing resilience of kelp beds along a latitudinal temperature gradient: potential implications for a warmer future. Ecology letters, 13(6), pp.685-694.
Becheler, R., Haverbeck, D., Clerc, C., Montecinos, G., Valero, M., Mansilla, A. and Faugeron, S., 2022. Variation in Thermal Tolerance of the Giant Kelp’s Gametophytes: Suitability of Habitat, Population Quality or Local Adaptation?. Frontiers in Marine Science, 9, p.802535.
Reed, D., Washburn, L., Rassweiler, A., Miller, R., Bell, T. and Harrer, S., 2016. Extreme warming challenges sentinel status of kelp forests as indicators of climate change. Nature Communications, 7(1), p.13757.
Rogers-Bennett, L. and Catton, C.A., 2019. Marine heat wave and multiple stressors tip bull kelp forest to sea urchin barrens. Scientific reports, 9(1), p.15050.